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Abstract 

Background: Genomic selection is an effective tool for animal and plant breeding, allowing effective individual 
selection without phenotypic records through the prediction of genomic breeding value (GBV). To date, genomic 
selection has focused on a single trait. However, actual breeding often targets multiple correlated traits, and, 
therefore, joint analysis taking into consideration the correlation between traits, which might result in more 
accurate GBV prediction than analyzing each trait separately, is suitable for multi-trait genomic selection. This would 
require an extension of the prediction model for single-trait GBV to multi-trait case. As the computational burden of 
multi-trait analysis is even higher than that of single-trait analysis, an effective computational method for 
constructing a multi-trait prediction model is also needed. 

Results: We described a Bayesian regression model incorporating variable selection for jointly predicting GBVs of 
multiple traits and devised both an MCMC iteration and variational approximation for Bayesian estimation of 
parameters in this multi-trait model. The proposed Bayesian procedures with MCMC iteration and variational 
approximation were referred to as MCBayes and varBayes, respectively. Using simulated datasets of SNP genotypes 
and phenotypes for three traits with high and low heritabilities, we compared the accuracy in predicting GBVs 
between multi-trait and single-trait analyses as well as between MCBayes and varBayes. The results showed that, 
compared to single-trait analysis, multi-trait analysis enabled much more accurate GBV prediction for low-heritability 
traits correlated with high-heritability traits, by utilizing the correlation structure between traits, while the prediction 
accuracy for uncorrelated low-heritability traits was comparable or less with multi-trait analysis in comparison with 
single-trait analysis depending on the setting for prior probability that a SNP has zero effect. Although the 
prediction accuracy with varBayes was generally lower than with MCBayes, the loss in accuracy was slight. The 
computational time was greatly reduced with varBayes. 

Conclusions: In genomic selection for multiple correlated traits, multi-trait analysis was more beneficial than 
single-trait analysis and varBayes was much advantageous over MCBayes in computational time, which would 
outweigh the loss of prediction accuracy caused by the approximation procedure, and is thus considered a 
practical method of choice. 
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Background 

A huge number of genome-wide polymorphisms have re- 
cently been elucidated in livestock and crops with the de- 
velopment of sequencing technologies. High-throughput 
genotyping systems, such as high-density SNP chips con- 
taining several tens or hundreds of thousands of genome- 
wide SNP markers and GBS (genotyping by sequence) [1], 
have become available to efficiently identify genotypes of 
individuals for a large number of SNPs at low cost. Conse- 
quently, genomic selection [2] is attracting attention as a 
new breeding technology utilizing the information of 
genome-wide dense SNP markers. In genomic selection, a 
model to predict genomic breeding value (GBV) based on 
genome-wide SNP genotype is firstly constructed using a 
training population consisting of individuals with both 
SNP genotypes and phenotypes of a trait and, subse- 
quently, using this model, the GBVs of a trait are predicted 
for individuals in the tested population, which are selec- 
tion candidates, based on their SNP genotypes, allowing 
effective individual selection to be performed without 
phenotypic records using the predicted GBV. Therefore, 
genomic selection requires the construction of prediction 
models being able to accurately relate genotypes of 
genome-wide SNPs to GBV. 

In the original study of genomic selection by Meuwissen 
et al. [2], two Bayesian methods were presented for con- 
struction of prediction models, which were referred to as 
BaeysA and BayesB and have been used for the studies of 
genomic selection in their original or modified forms [3,4]. 
The model of BayesA is equivalent to the Bayesian shrink- 
age regression (BSR) model [5], where all SNPs are 
included in the prediction model as covariates and the esti- 
mates of SNP effects are shrunk by assuming a normal dis- 
tribution with mean 0 and SNP-specific variance as prior 
distributions of SNP effects, resulting in negligible estimates 
for the effects of many SNPs irrelevant to phenotype, but 
significantly large estimates for the effects of SNPs contrib- 
uting to phenotype. On the other hand, the BayesB method 
is regarded as a variant of Bayesian stochastic search vari- 
able selection (BSSVS) [6], where the prior probability, n, 
that a SNP has zero effect and is removed from the model 
was considered in the model fitting to obtain the best 
model explaining the phenotypes with a small number of 
SNP effects. A normal distribution with mean 0 and SNP- 
specific variance is assumed for the prior distribution of 
SNP effects in BayesB, as in BayesA, if SNPs are included 
in the model with probability 1-n and, otherwise, both the 
mean and variance are fixed at 0 with probability n. 

In BayesA and BayesB, an additional hierarchical 
structure is induced for this SNP-specific variance, 
where an inverted chi-square distribution with degree of 
freedom v and scale parameter S is adopted as the prior 
distribution of the variance. However, only the informa- 
tion of the relevant single SNP can be used for the 



posterior inference of this SNP-specific variance regardless 
of the number of genotypes or phenotypes. Due to the in- 
sufficient Bayesian learning for the variance, the degree of 
shrinkage for the estimates of SNP effects is largely influ- 
enced by the prior setting for v and S in BayesA and 
BayesB [7]. In BaeysB, the sparseness of the model and the 
prediction accuracy are also greatly affected by the prior 
probability, n, that a SNP has zero effect, which is treated 
as a known fixed value. To overcome these drawbacks of 
BayesA and BayesB, the modified methods such as BayesC 
and BayesD were proposed. In BayesC, a single variance 
common to all SNPs is adopted for the prior distribution 
of SNP effects while BayesD allows S to be inferred from 
the data, given v. Furthermore, the step for inferring n, 
which is given as a fixed value in BayesB, can be incorpo- 
rated in the procedure of BayesC and BayesD and the cor- 
responding methods are termed BayesC/r and BayesD7r, 
respectively [4]. 

These Bayesian methods were mainly developed for gen- 
omic selection of a single trait. However, actual breeding of 
animals and plants often aims to simultaneously improve 
multiple correlated traits. Therefore, joint prediction of 
GBVs for multiple traits, taking into consideration the cor- 
relation structure between traits, is suitable for multi-trait 
genomic selection, which requires the extension of existing 
methods for single-trait GBV prediction to multi-trait case. 
In QTL mapping methods that use similar models to those 
of genomic selection, some researchers have developed 
multi-trait models [8-10]. Xu et al. [8] extended the BSR 
model to a multi-trait case introducing the correlation 
structure between traits in estimation of QTL effects and 
other non-genetic effects. Banerjee et al. [9] employed a 
Bayesian composite model space approach [11,12] for 
multi-trait QTL mapping, where a variable indicating inclu- 
sion or exclusion of each QTL was incorporated for con- 
structing a model with the smallest possible number of 
QTLs, a similar approach to BSSVS adopted in BayesB, and 
each trait was allowed to have a different model by assum- 
ing the trait-specific effects for QTLs and non-genetic fac- 
tors that were assumed to be uncorrelated between traits. 
Meuwissen and Goddard [10] proposed a multi-trait BSSVS 
model for QTL mapping. These methods can be applied 
for prediction of multi-trait GBVs in genomic selection. 
Calus and Veerkamp [13] applied the method proposed in 
[10] with some modification for multi-trait GBV prediction 
to compare the accuracy between single-trait prediction 
and multi-trait prediction in genomic selection. 

The computational procedure with MCMC iteration is 
generally used for the Bayesian methods to estimate 
parameters in the models, which become complicated 
models in genomic selection, including a huge number 
of SNPs as covariates and SNP effects that are estimated 
as their regression coefficients. The computational bur- 
den of MCMC-based Bayesian methods, which requires 
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a long time until convergence when estimating many 
parameters is huge even in single-trait GBV prediction 
and would be further increased in the case of multiple 
traits, thus hindering the MCMC procedure depending 
on the number of traits to be jointly analyzed. Therefore, 
it would be necessary to devise a solution for reducing 
the computational burden of Bayesian methods for 
multi-trait GBV prediction. So far, some non-MCMC 
computational procedures for Bayesian methods have 
been proposed in QTL mapping and genome-wide asso- 
ciation study, including EM- algorithm [14] and vari- 
ational approximation [15,16]. The EM-algorithm was 
also applied to single-trait GBV prediction, successfully 
reducing the computational time [17,18]. 

In this paper, we propose Bayesian methods for multi- 
trait GBV prediction, in which BSR models allowing 
variable selection are developed and MCMC procedures 
for estimating model parameters including SNP effects 
are described as well as a computationally cost-effective 
non-MCMC method using variational approximation as 
an alternative computational procedure. Hereafter, the 
Bayesian methods based on MCMC iteration and vari- 
ational approximation are referred to as MCBayes and 
varBayes, respectively. The multi-trait Bayesian models 
described here include a Bayesian shrinkage regression 
(BSR) models that are equivalent to those adopted by 
BayesA and BayesD when variable selection is not con- 
ducted, and the models of BSSVS that are regarded as 
slightly modified versions of BayesB and BayesD7r meth- 
ods when a step of variable selection step is incorpo- 
rated. We develop computationally-effective variational 
approximation procedures to construct the posterior dis- 
tributions of parameters of these Bayesian models. 

Using simulated datasets consisting of genotypes of 
genome-wide dense SNPs and phenotypes of three corre- 
lated traits with high and low heritabilities, we investigated 
the differences in prediction accuracy for each trait be- 
tween multi- trait analysis, where GBVs of three traits were 
simultaneously predicted taking the correlation structure 
between traits into consideration, and single-trait analysis, 
where each trait was separately predicted for GBV. We 
also evaluated the prediction accuracy of the varBayes 
methods in comparison with MCBayes. Moreover, we 
investigated the performance of multi-trait analysis in 
simulated data including missing phenotypes. 

Methods 

In this section, we describe Bayesian models for multi- 
trait GBV prediction and computational procedures for 
Bayesian estimation of the model parameters, including 
construction of posterior distributions of the parameters 
using MCMC iteration and variational approximation. 
Here, we consider the statistical models for BSR and 
BSSVS, which are shown to be equivalent to BayesD and 



similar to BayesD7T, respectively. In these Bayesian mod- 
els concerned, BSR is regarded as a special version of 
BSSVS by setting n =0, where n is a prior probability 
that a SNP does not contribute to traits; thus, we 
present the multi-trait GBV prediction models in a uni- 
fied fashion in terms of BSSVS. 

We assume that the number of SNPs genotyped is N 
and a training dataset including n individuals with 
phenotypic records of multiple traits, where the number 
of traits is denoted by T, and SNP genotypes is available 
for estimating parameters in the prediction model. We 
also assume that a test dataset consists of individuals 
with only SNP genotypes, for each of which GBV is pre- 
dicted. We denote two alleles at each SNP by 0 and 1 
and three genotypes by 'Q_Q', '0_1', and T_T. 

Models for Bayesian stochastic search variable selection 
in multi-trait genomic selection 

We propose the following Bayesian multi-locus linear 
model for the phenotypes of T traits of the ith individ- 
ual, which are denoted as a 7x1 vector y, (/ = 1,2,. . .,«), 
as a BSSVS model for multi-trait GBV prediction: 



Yi = X,b + Yi u n%i + e < 



(1) 



where b is a/x 1 vector of non-genetic effects including 
interceptions of the model with X, being a T xf design 
matrix linking b to the ith individual; w a is a variable in- 
dicating the genotype of the ith individual at the /th SNP 
taking a value of -1, 0 or 1 corresponding to the geno- 
types, '0_0! '0_1J or '11! respectively; g/ is a Tx 1 vector 
of the effects of the /th SNP on the phenotypes of T 
traits; Yl is a variable indicating the inclusion of the /th 
SNP in the model with 1 or 0 depending on whether or 
not the /th SNP is included in the model; and e, is a T x 
1 vector of residual errors following a T-variate normal 
distribution N(0, E e ) with 0 being a T x 1 vector of zeros 
and Yj e being a TxTcovariance matrix of e,. 

Within the Bayesian framework, prior distributions are 
assigned to the parameters of the model (1). We assume 
that the priors of the elements of b are the improper 
uniform distribution over the possible values. The prior 
probabilities that yi =1 an d Yi =0 are presented as 1- n 
and Ji, respectively, considering the prior probability that 
a SNP does not contribute to the trait and is excluded 
from the model is n. The prior distribution of g; given y c 
has the form 



&i\Yf 




(2) 



where ~L g i is a T x T matrix that is the variance and co- 
variance matrix of g/ given yi =1 and d (0) is the delta 
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function that concentrates a total mass at zero for all T 
elements of g/. We induce a hierarchical structure for ~L g i 
by assigning an inverse Wishart distribution with degree 
of freedom v and scale parameter S, denoted by IWj{v,S), 
to the prior distribution of Y. g i : 



£ s /~iWr(v,S)oc|s| v/2 |^ 



|-(v+T+l)/2 



exp{--tr(S^r 1 )} 
(3) 



Although the SNP effect g; is zero and irrelevant to Y. g[ 
when Yi =0 as shown in (2), we assume that Y. g i is a 
priori distributed as given in (3) regardless of yi- We 
treat v as a given fixed value but infer S from the data 
within the Bayesian framework. For simplicity, we as- 
sume here that S is a diagonal matrix with the /<th diag- 
onal element being s k , that is S = diag^) (/c= 1,2,. . .,T), 
and we adopt the improper uniform distribution over 
positive values for a prior distribution of s k . The residual 
variance and covariance matrix S e is assumed to have a 
uniform distribution over the positive definite matrices 
as a prior distribution. 

Denoting these parameters of the Bayesian model col- 
lectively by 8, the prior distribution of 0 by p(&) and 
observed phenotypes y,- and SNP genotypes Uu (i = 1,2,. . ., 
n;l = 1,2,. . .,N) by Y and U, respectively, we can write the 
posterior distribution of 0, g(B\ v,Y, U), as 



^|v,Y,U)oc|j: e |-"/ 2 e X p(-^ yi *% 
z i=i 



1. 



'f[[{ (1-^)1^1 



-1/2 



expt-ig/iV^A^O)} 1 - 



<ni s i v/a i^ 



i-(v+r+i)/2 



exp( 



-trS£„ 



(4) 



from (1), (2) and (3), where y* is a residual given by 

N 

Yi* = Yi ~ x < b " Yi u ^i- 
l=i 

Given that S is fixed, posterior distribution (4) is 
equivalent to that of BayesA extended to multi-trait case 
when 77 = 0 leading to ji =1 f° r all SNPs. When 7r>0, 
the selection of SNPs to be included in the model is per- 
formed as in BayesB. However, it should be noted that 
posterior distribution (4) is not equivalent to that of the 
multi-trait version of BayesB, which supposes that is 
a zero matrix as well as gi when yi =0, while, in (4), the 
prior distribution of Y. g i is assumed to be IWj(v, S) re- 
gardless of yi- The form of posterior distribution (4) 
allows Gibbs sampling in the MCMC estimation for all 
parameters including gi and ~L g i. The full conditional 
posterior distributions of parameters used in MCBayes 
are described in Additional file 1. 



Variational approximation procedure for multi-trait 
Bayesian model 

We adopted variational approximation as an alternative 
to MCMC iteration for constructing marginal posterior 
distributions of parameters based on joint posterior dis- 
tribution (4). In the variational approximation proced- 
ure, the joint posterior distribution is approximated by 
the product of functions for subsets of parameters with 
lower dimension. Briefly, we assume that g(B 1 v, Y, U) is 
approximated by a function of 6, q{B), which is factor- 
ized as q(B) = 171(81) #2(62)- • ?jc(9jc). where 81, 8 2) . . . ,B K 
are mutually exclusive subsets of 6 such that U k = 18^ = 8 
and qk (B k ) (k = 1,2,. . .,K) may generally depend on v,Y 
and U although this dependence is omitted here for sim- 
plicity. This approximating marginal posterior distribu- 
tion, qiJJdk): is referred to as the variational posterior of 
6/r [19]. The form of qk{Qk) is determined such that the 
Kullback-Leibler divergence between g(8|v,Y, U) and q 
(8 1 tfi(Bi) #2(82). • -fete*), D(q\\g), is minimized [20], 
where D(q \ \g) is defined as 



D(q\\g)= / 4(8) log 



q{%) 



'£(8|v,Y,U) 



dB 



It can be shown that ^(8^) is expressed as 
q k (Q k ) = Cexp{£_8*[log£(8|v,Y,U)]} (5) 

where C is a normalized constant and E^ ej t[-] indicates 
the expectations of parameters other than 8^ that are 
calculated with respect to every other parameter's vari- 
ational posteriors except ^(8^.) [21]. 

In the varBayes method that applies the variational ap- 
proximation procedure to the multi-trait Bayesian model 
considered here, g(8|v,Y, U) is assumed to be approxi- 
mated by a factorized function ^(6) that is written as 

TV 

q(d) = q(b)q(E e )l[ {?(r„ft)?(^)}?(S) 

1=1 

where we denote all of the variational posteriors for the 
different parameters in the right-hand side by q(.) for 
simplicity with the understanding that q{b), q (S e ) and 
so on take different forms depending on the parameters. 
The forms of these variational posteriors are derived 
from (5) in a manner similar to that used for derivation 
of the full conditional posteriors for the parameters in 
Gibbs sampling and the computational details are given 
in Additional file 2. In the following, we will give 
the variational posteriors for parameters, yi> Si> ^gi 
(1=1,2,. . .,7V), b, Z e and S. 

The variational posterior for b, q(b),is a multivariate 
normal density with mean 
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£(b) = {^X/£(Z e - 1 )X i }- 1 J;x/£(5: e - 1 ) 

1=1 n f'=l 

iy* -£>*/£ (r/&)} (6) 



/-l 



and variance-covariance matrix 

y(b) = {^x/£(z e - 1 )x,}- 



where £(.) indicates the expectation calculated with re- 
spect to the variational posteriors of relevant parameters, 
while q (E e ) is IW T (n-T-1, S e ) with 



N 

-i 



s * = E{y.- X < £ ( b ) - E 

i=i /=i 

TV 

{ yi -Xi£(b)-^ Ma £( r/ g,)}' (7) 

from which we obtain 

Efa- 1 ) = (n - T - lJST 1 (8) 

The variational posterior of 5^ is represented by IW r 
(v^, S^), where 

iy = v + E( Yl ) 

and 

= £(S)+£(/,ftg/) 

The expectation of Y. g i ' with respect to this posterior 
distribution is 

Efigt *) = v giSgi 1 

= {v + £( r/ )}{£(S)+£( r/ g,g/)} (9) 

For Yi an d gz> we consider joint distribution for vari- 
ational posterior q(yi, El) that is expressed as 

q( Yl , ft) oc[(l - ^V/^Sr 1 ) | 1/2 exp{^ g 'V,,- 1 g } 



«p{gi\gi,Vgi)) r 'x[nS(0). 



(10) 

where </>(g/| g/,Vg/) is a density function of a multivari- 
ate normal distribution /V(g ; , V^) with mean 

w n 

ft = ^(zr 1 ) +^^ 2 £(£ e " 1 )}" 1 £(£ e " 1 )^« i / 
i=i i=i 

{ yi -X,E(b)-53«^E(y^)} (11) 

/*' 

and variance-covariance matrix 



v^ = {£(i:r 1 ) + E M « 2£ (^" 1 )}" 1 



(12) 



The marginal posterior distribution of 77 is a binomial 
distribution with probabilities of Y i =1 an d 0 given by 

?(rz = i) = E(r,) 

(i - ^iv^i^iECzr 1 )! 1 ^ «p{i g'^r 1 &} 



and 



(i - w )|v„HE(Zri)|v»«p{i g^r 1 g ; } + n 

(13) 

?fo = o) = i - = i) 



(l-7r)|V^|V2|£(s r iMV2 



exp{| g^Vji : ft} +7T 
The conditional distribution of g/ given yj is given by 



6(0) 



(r/ = o) 



Therefore, we obtain 

= = !) £ (gz|Xz = !) + ?(y« = 0)E( gl \ Yl = o) 
= £ (X/)g/ 



and 



E{YlSlSl) = E {Yi) ( 8/ 8/ + v ? 



(14) 



(15) 



From (4), the variational posterior of S is a £-variate 
Wishart distribution with a scale matrix Z s and degree 
of freedom v s , W> s , E s ), where I, = i^ 1 )} - 1 
and v s = A/V + £+ 1, and the expectation of S is expressed 
as 

£(S) = v s E s 

= (wv+r + i){^f = £(s,r 1 )}-i (i6) 

As outlined above, a well-known probability distribu- 
tion, such as normal, inverse Wishart and so on, is 
assigned to the variational posterior of each parameter, 
which is characterized by the expectations of the func- 
tions of other parameters, taken with respect to their 
variational posteriors. The relationships between these 
expectations are given by (6), (8), (9), (13), (14), (15) and 
(16), from which the expectations can be calculated with 
numerical iterations to obtain the variational posteriors. 

Moreover, the prior probability for a SNP to have zero 
effect, Ji, can be treated as a variable parameter to be in- 
ferred from the data as in BayesC/r and BayesD7r when 
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0 < Ji < 1. Here, we assume a uniform prior on 0 < n < 1 
to obtain a Beta distribution, B(a,b), as a variational poster- 
ior from (4), where a = N- £ ?= iE{yD + 1 and b = £ f = X E 
(yi) + 1, with the expectation 

E(n) = a/(a + b) 

= { 7V -EL £ W + 1 }/( Ar + 2 ) ( iy ) 

Accordingly, the variational posteriors of the other 
parameters are modified with n substituted by £(77) 
when ji is involved in the Bayesian inference. 

Treatment of missing phenotypes 

In the phenotypic records for multiple traits, it is com- 
mon for trait values to be partially missing in some indi- 
viduals. Missing phenotypes of a trait in individuals can 
be inferred with the observed phenotypes of other traits 
in the same individuals. The step for inferring missing 
phenotypes can be implemented in the MCBayes and 
varBayes procedures as described below. 

When there are missing phenotypes in an individual, 
the residual vector of the individual, e in model (1), is 
partitioned into components e Q and e m corresponding to 
observed and missing phenotypes, respectively. Follow- 
ing [22], e m can be sampled from the following normal 
distribution, 

where 2,^, E mo and S 00 indicate the partition of the re- 
sidual variance-covariance matrix ~L e corresponding to e m 
and e 0 . Accordingly, e m is drawn with Gibbs sampling in 
MCBayes while it is obtained as E(Z mo )E(E 00 )" 1 e 0 * in var- 
Bayes, where, denoting the component corresponding to 
the missing traits by subscript 'o,' e Q * is written as 

N 

e 0 * = y 0 - X 0 £(b) - u u E (ViSoi) 

1=1 

and the expectations can be calculated with the variational 
posteriors of E e , b and gi . Missing phenotypes are inferred 
as the sum of the estimates of b, g ; and e m , which are used 
for the construction of the prediction model. 

Simulation experiments 

We simulated datasets to evaluate the accuracy of the 
predicted GBVs using the proposed Bayesian methods, 
MCBayes and varBayes, for multiple traits. In generating 
the datasets, three traits, denoted as A, B and C, were 
considered. The simulation of population and genome 
was carried out following [3] where a single trait was 
treated, while multiple traits were generated through a 
modified approach. The simulated population, genotypes 
and phenotypes are described in the following. 



Populations with an effective population size of 100 
were maintained by random mating for 1000 generations 
to attain mutation drift balance and LD between SNPs 
and QTLs. In generation 1001 and 1002, the population 
size was increased to 1000. The population in the 1001st 
generation was treated as a training population, where 
the phenotypes of three traits and SNP genotypes of the 
individuals were simulated and analyzed to estimate the 
SNP effects in the model. The phenotype of each trait 
for each individual in the 1001st generation was given as 
the sum of QTL effects over the polymorphic QTLs and 
environmental effects, which were sampled as described 
later. For simplicity, no other fixed effects were assumed. 
The population in the 1002nd generation was used as a 
test population, where the individuals were only geno- 
typed for SNP markers without phenotypic records and 
GBVs of three traits were predicted for each individual 
using a model with SNP effects estimated based on the 
training population in the 1001st generation. The true 
breeding value (TBV) of the individual in the 1002nd 
generation was also simulated as the sum of QTL effects 
corresponding to the QTL genotype for each trait and 
used for evaluating the accuracy of predicted GBV, but 
was regarded as unknown and unavailable in the estima- 
tion of SNP effects in the models. Accuracy was mea- 
sured based on the correlation between the TBV and 
predicted GBV, >"tbv,pGbv. for each trait and regression of 
TBV on predicted GBV, &tbv, p gbv> was also obtained for 
assessing the bias of the prediction. 

The genome was assumed to consist of 10 chromo- 
somes each lOOcM in length. Two scenarios were consid- 
ered for the number of available SNP markers and the 
datasets under these two scenarios were denoted as Data I 
and Data II. In Data I, 101 marker loci were located every 
lcM on each chromosome for a total of 1010 markers on 
a genome. In Data II, 1010 equidistant marker loci were 
located on each chromosome for a total of 10100 markers. 
We assumed that 100 equidistant QTLs were located on 
each chromosome such that a QTL was in the middle be- 
tween two marker loci in both Data I and Data II. There- 
fore, there were a total of 1000 QTLs located on a whole 
genome. The mutation rates assumed per locus per mei- 
osis were 2.5 x 10" 3 and 5.0 x 10' 5 for the marker locus 
and QTL, respectively. At least one mutation occurred in 
the most of the marker loci with a high mutation rate dur- 
ing the simulated generations. In the marker loci experien- 
cing more than one mutation, the mutation remaining at 
the highest minor allele frequency (MAF) was regarded as 
visible, whereas the others were ignored, which resulted in 
the marker loci having two alleles similar to SNP markers. 
Although the mutation rate for QTL was assumed 2.5 x 
10' 5 in the simulation for a single trait conducted in [3], 
we here doubled it for generating TBVs of multiple traits 
for the reason as described below. 
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The polymorphic QTLs at which mutation occurred 
were used to simulate the three traits, A, B and C, the her- 
itabilities of which, denoted by h\, h\ and h c , respectively, 
were assumed to be h\ = 0.8, /?b = 0.1 and he = 0.1. We 
divided all polymorphic QTLs into four groups, Groupl, 
Group2, Group3 and Group4, according to the causal re- 
lationship with traits, where Groupl was a group of QTLs 
affecting both traits A and B, Group2 was a group of 
QTLs affecting only trait A, Group3 was a group of QTLs 
affecting only trait B and Group4 was a group of QTLs 
affecting only trait C. Therefore, it was assumed that 
QTLs of Groupl had pleiotropic effects on traits A and B 
while QTLs of other groups affected only a single trait. 
Each polymorphic QTL was randomly assigned to one of 
these groups, that is, to Groupl, Group2, Group3 and 
Group4 with respective probabilities 0.4, 0.1, 0.1 and 0.4. 
Therefore, 80% of the QTL loci affecting trait A were 
shared by trait B on average and whereas all the QTL loci 
affecting trait C influence neither traits A nor B. In this 
setting of multi-trait case, the mutation rate of QTL was 
increased to 5.0 x 10' 5 from 2.5 x 10" 5 , the mutation rate 
adopted in [3] for generating a single trait, to retain the 
number of QTL per trait. 

The pleiotropic effects of each QTL in Groupl were 
assumed to be correlated between traits A and B with 
correlation coefficient of 0.9. Consequently, genetic cor- 
relations between traits A and B, between A and C and 
between B and C, denoted as /3gab< Pgac and pese, re- 
spectively, were p GAB = 0.72 and /3 G ac = £gbc = 0 on 
average although the values of correlation coefficients 
were somewhat fluctuated in each data generation. This 
setting of traits was adopted to investigate how the pre- 
diction accuracy was increased for a low-heritability trait 
(trait B) by simultaneous analysis for multiple traits in- 
cluding a correlated high-heritability trait (trait A). 

The effects of QTL alleles were sampled from gamma 
distributions independently for each QTL. Pleiotropic 
effects of QTL alleles in Groupl were determined for 
traits A and B by generating two correlated gamma ran- 
dom variables, x and y, with the correlation coefficient of 
0.9, and assigning them with positive or negative values 
with equal probabilities. The effects of QTL alleles in the 
other groups, which affected each single trait only, was 
determined by sampling a gamma random variable z and 
by similarly assigning it with the positive or negative sign. 
We generated these three random variables x, y and z such 
that their marginal distributions were the same gamma 
distribution with scale parameter 0.4 and shape parameter 
1.66, Gamma(0.4,1.66). For obtaining correlated variables 
x and y, we generated three independent gamma variables, 
X\ } X2 and #3, which were sampled from Gamma 
(0.36,1.66), Gamma(0.04,1.66) and Gamma(0.04,1.66), re- 
spectively, and determined the values of x and y as x = X\ 
+ x 2 and y = Xi+ X3. It can be shown that x and y had a 



correlation coefficient of 0.9 and the same marginal distri- 
bution Gamma(0.4,1.66) [23]. 

The environmental correlation coefficients between 
three traits were denoted by /3eab» Peac and /?ebc> re- 
spectively, and assumed as /?eab = 0.1, Peac = 0.2 and 
Pebc = 0-3- The environmental effects were sampled from 
a trivariate normal distribution with a mean vector 0 
and a variance-covariance matrix R E , where 

(°EA 2 ^eAbOEAOEB £eAC°EAOeC \ 

/>eab°ea0eb Ceb 2 /'ebc^ebO'ec I 
ftAC^EAffEC ^EBcOEBffEC ^ / 

with Oea, Oeb and oic indicating environmental variances 
of three traits. Environmental variance of trait A was 
given by oea = (1/^a-1)oga with its heritability h\ and 
genetic variance Oga> which was variance of TBV of trait 
A, and those of traits B and C were obtained similarly. 
The environmental effects were added to TBVs which 
were given by the sum of QTL effects to determine 
phenotypic values of three traits for individuals in the 
1001st generation. 

In Data I, 100 replicated datasets were simulated while 
Data II consisted of 20 replicated datasets due to the larger 
number of SNPs. Each of replicated datasets included 
records of phenotypes of three traits and genotypes of 
SNPs for the training population (1001st generation) and 
only SNP genotypes for the test population (1002nd gen- 
eration). To simulate the situation of missing phenotypes, 
we generated additional datasets by deleting the pheno- 
typic records of some traits for some individuals in the 
100 replicated training datasets of Data I. These 100 repli- 
cated datasets were referred to as Data III, where the 
phenotypic records of traits A, B and C were respectively 
deleted for individuals of i = 801-1000, individuals of 
i = 1-500 and individuals of i = 201-700 in 1000 individuals 
{i = 1-1000) of the 1001st generation of Data I. Therefore, 
in Data III, the prediction model for GBVs was con- 
structed with a training dataset consisting of the pheno- 
types of 800, 500 and 500 individuals for traits A, B and C, 
respectively, where only 100 individuals (z = 701-800) had 
phenotypic records of all three traits, and 1000 individuals 
in the 1002nd generation were predicted for GBV based 
on the genotypes of 1010 markers. The setting of non- 
zero environmental correlations, i.e., Peab = 0.1, Peac = 0-2 
and Pebc = 0.3, was adopted here to assess the benefit from 
the estimation process of missing phenotypes implemented 
in multi-trait analyses for the prediction accuracy, where 
the information of environmental covariance between 
observed and missing phenotypes was utilized as in (18). 

Each replicated dataset in Data I, Data II and Data III 
was analyzed using the proposed methods for multiple 
traits, MCBayes and varBayes, to construct the GBV 
prediction model in thelOOlst generation and investigate 
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^tbv.pGbv and brsv, p gbv for each trait in the 1002nd 
generation. For a comparison with conventional single- 
trait GBV prediction, the same datasets were also ana- 
lyzed with the single-trait setting of MCBayes where 
three traits were separately treated without considering 
the correlation structure between traits. 

We conducted cross-validation as well to evaluate the 
prediction performance within a population in the 1001st 
generation without the 1002nd generation since the techni- 
ques of cross-validation have commonly been used for 
evaluation of the accuracy in the studies of genomic selec- 
tion with the actual datasets of animals [24-26] and plants 
[27-29], where the prediction accuracy of the model for the 
unobserved future samples was of concern. We applied 10- 
fold cross-validation. In brief, we randomly split a popula- 
tion of 1000 individuals in the 1001st generation into ten 
subpopulations with each size 100. By using a single subpo- 
pulation used as a test set and the remaining nine subpopu- 
lations as a training set to construct prediction model, the 
prediction accuracy and bias were assessed for the test set. 
This process was repeated ten times until all single subpo- 
pulations were used as test sets exactly once. We averaged 
the prediction accuracy and bias over ten repetitions to 
evaluate the prediction performance in each dataset. Be- 
cause of much computational burden with MCMC analysis 
for repeated model constructions in cross-validation, evalu- 
ation with 10-fold cross-validation was carried out for Data 
I only, where mean of the prediction accuracies and biases 
in 100 datasets were obtained as r TB v, P GBV and £>tbv, p gbv 
for each trait and for each prediction method. In the 
process of cross-validation, we also investigated the correl- 
ation between predicted GBV and the phenotypic value, 
r ypGBV , in place of TBV, and regression of predicted GBV 
on phenotypic value, b Yi p gbv> for each trait. 



Two settings for the prior probability that a SNP has 
zero effect, n, were adopted in model construction, i.e., n 
was fixed at 0 or n was variable over the range 0 < n <1 
and inferred from the data. The values of hyperparameter 
v, degree of freedom of prior inverse Wishart distribution 
of Eg/, were set to 5.0 and 3.2 corresponding to n =0 and 
0 < 77 <1 after preliminary analyses to evaluate the effects 
of the values of v on prediction accuracy over 3 < v <6 
although the accuracy was little affected by the values of v. 

In the MCMC iteration of MCBayes, we repeated 11000 
cycles including a burn-in period of the first 1000 cycles. 
The values of parameters were sampled every 10 cycles to 
obtain the posterior means that were used to determine a 
prediction model for each generated dataset. In the method 

of varBayes, we adopted the criterion | 8 — 8| 2 /| 0 | 2 < 10~ 8 
for convergence in numerical iteration for computing the 

expectations of parameters, where 8 and 8 are the current 
and previous value of the expectations of parameters. When 
this criterion was satisfied, the computational procedure 
with variational approximation was regarded as converged. 

Results 

We evaluated the accuracy and bias of the predicted 
GBVs, r TBV ,pGBv and &tbv, p gbv« obtained by the pro- 
posed methods for genomic selection of multiple traits 
including MCBayes and varBayes in 100 simulated data- 
sets of Data I and Data III and in 20 datasets of Data II 
in comparison with the prediction accuracy with single- 
trait MCBayes analysis. The means of accuracies and 
biases of prediction evaluated using a test population in 
the 1002nd generation are summarized in Table 1, 
Table 2 and Table 3 for Data I, Data II and Data III, re- 
spectively. In all of datasets, MCBayes provided higher 



Table 1 Accuracy and bias of predicted GBVs in Data I 



Method 






Trait A 


Trait B 


Trait C 


MCBayes 


77 = 0 


'TBV.pGBV 


0.788 ±0.051 


0.581 ±0.103 


0.453 ± 0.090 






^TBV.pGBV 


0.994 ±0.038 


1 .048 ± 0.264 


1.00 ±0.370 




0<77 <1 


'TBV.pGBV 


0.753 ± 0.060 


0.580 ±0.1 17 


0.364 ±0.1 37 






^TBV.pGBV 


1 .070 ± 0.064 


1.1 49 ±0.340 


1.01 6 ±0.364 


varBayes 


77=0 


'TBV.pGBV 


0.754 ±0.061 


0.570 ±0.1 13 


0.383 ±0.1 17 






^BV.pGBV 


1.054 ±0.051 


0.994 ± 0.233 


0.899 ± 0.247 




0<77 <1 


'TBV.pGBV 


0.71 6 ±0.070 


0.548 ±0.1 22 


0.347 ±0.1 31 






^TBV.pGBV 


0.894 ±0.054 


0.834 ±0.1 86 


0.636 ± 0.202 


single-trait 


77=0 


'TBV.pGBV 


0.783 ±0.051 


0.469 ± 0.083 


0.455 ± 0.076 


(MCBayes) 




^TBV.pGBV 


0.978 ± 0.037 


1.020 ±0.301 


0.970 ± 0.259 




0<77 <1 


'TBV.pGBV 


0.778 ±0.050 


0.491 ±0.114 


0.483 ±0.101 






^BV.pGBV 


1 .089 ± 0.054 


1.1 10 ±0.634 


1.061 ±0.338 


Averages and standard errors based on 
probability that a SNP has zero effect, n, 


100 replicates of simulated data 
we considered two settings, in 


are listed for prediction accuracy, r pgbv.tbv, an d bias, bpGBv/rev. of each trait. For the prior 
which n was fixed at 0, meaning the inclusion of all SNPs in the model, and n was varied 



over 0 < jt <1 and inferred from the data. 
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Table 2 Accuracy and bias of predicted GBVs in Data II 


Method 






Trait A 


Trait B 


Trait C 


MCBayes 


77 = 0 


r TBV,pGBV 


0.902 ± 0.032 


0.706 ±0.1 03 


0.519 + 0.097 






frlBV.pGBV 


0.998 ± 0.034 


0.902 ±0.1 11 


0.796 ±0.1 79 




0<77 <1 


tBV.pGBV 


0.868 ± 0.047 


0.731 ±0.120 


0.401 ±0.182 






^TBV.pGBV 


1 .092 ± 0.093 


1.1 89 ±0.1 99 


1.1 98 ±0.553 


varBayes 


77 = 0 


'YBV.pGBV 


0.859 ± 0.049 


0.656 ±0.1 10 


0.438 ± 0.074 






bTBV.pGBV 


1 .059 ± 0.065 


0.799 ±0.1 05 


0.724 ±0.1 11 




0 < 77 <1 


'YBV.pGBV 


0.838 ±0.061 


0.678 ±0.140 


0.330 ±0.1 57 






faTBV.pGBV 


0.983 ± 0.034 


0.851 ±0.138 


0.562 ±0.1 55 


single-trait 


77=0 


'YBV.pGBV 


0.884 ±0.039 


0.485 ± 0.086 


0.493 ± 0.089 


(MCBayes) 




^BV.pGBV 


0.974 ±0.035 


0.766 ±0.1 13 


0.766 ±0.1 13 




0 < 77 <1 


(YBV.pGBV 


0.843 ± 0.044 


0.597 ±0.1 20 


0.601 ±0.109 






bTBV.pGBV 


1.562 ±0.261 


1.787 ±0.431 


1.832 ±0.565 



Averages and standard errors based on 20 replicates of simulated data are listed for prediction accuracy, r pG Bv.7Bv, and bias, bpGBv.TBv, of each trait. For the settings 
of the prior probability that a SNP has zero effect, tt, see Table 1 . 



prediction accuracy than varBayes in multi-trait analysis 
(Tables 1, 2 and 3). In Data I, r TBV)pGBV obtained with 
multi-trait MCBayes analysis was 0.788 (0.051), 0.581 
(0.103) and 0.453 (0.090) for traits A, B and C, respect- 
ively, when ji was fixed at 0, where standard errors were 
given in parenthesis here and hereafter, while that was 
0.753 (0.060), 0.580 (0.117) and 0.364 (0.137) when n 
was varied and inferred. In Data II with the number of 
SNPs increased to 10100, J" TBV , p gbv with multi-trait 
MCBayes analysis was higher at 0.902 (0.032), 0.706 
(0.103) and 0.519 (0.097) for traits A, B and C, respect- 
ively, when n was fixed at 0 while that was 0.868 (0.047), 
0.731 (0.120) and 0.401 (0.182) when n was inferred. 
With multi-trait varBayes analysis, r TBV>P GBv was 0.754 
(0.061), 0.570 (0.113) and 0.383 (0.117) for traits A, B 
and C with n =0 and that was 0.716 (0.070), 0.548 



(0.122) and 0.347 (0.131) with n variable in Data I while 
that was 0.859 (0.049) , 0.656 (0.110) and 0.438 (0.074) 
with ji =0 and 0.838 (0.061) , 0.678 (0.140) and 0.330 
(0.157) with 77 variable in Data II. The prediction with 
multi-trait MCBayes was almost unbiased in Data I, 
where &tbv, p gbv was near 1, but was biased in Data II 
for traits B and C (Tables 1 and 2). The analysis with 
varBayes showed greater bias than with MCBayes. 

In single-trait analysis where MCBayes was applied for a 
single-trait model, r TB v, P GBV was comparable with that of 
multi-trait analysis for high-heritability trait A while it was 
significantly decreased for low-heritability trait B which 
was highly correlated with trait A (Table 1 and Table 2). 
For trait C which had same heritability as trait B but was 
genetically independent of trait A (/3gac = 0.0), multi-trait 
MCBayes analysis gave the accuracy comparable to single- 



Table 3 Accuracy and bias of predicted GBVs in Data III 



Method 






Trait A 


Trait B 


Trait C 


MCBayes 


77=0 


'YBV.pGBV 


0.766 ±0.058 


0.500 ±0.1 27 


0.322 ±0.082 






^TBV.pGBV 


0.977 ±0.048 


0.998 ± 0.356 


0.967 ± 0.773 




0 < 77 <1 


fTBV.pGBV 


0.723 ±0.069 


0.503 ±0.141 


0.202 ±0.1 19 






£>TBV,pGBV 


1.065 ±0.076 


1.1 95 ±0.530 


0.799 ±0.523 


varBayes 


77=0 


r TBV,pGBV 


0.726 ± 0.072 


0.447 ± 0.131 


0.261 ±0.134 






^-TBV.pGBV 


0.984 ±0.052 


0.582 ±0.241 


0.383 ±0.1 81 




0 < 77 <1 


''TBV.pGBV 


0.679 ±0.081 


0.387 ±0.1 15 


0.228 ±0.1 12 






^TBV.pGBV 


0.840 ± 0.068 


0.389 ±0.1 32 


0.240 ±0.1 10 


single-trait 


77=0 


r TBV,pGBV 


0.760 ±0.058 


0.345 ± 0.070 


0.336 ±0.068 


(MCBayes) 




OTBV,pGBV 


0.965 ± 0.047 


0.931 ±0.368 


0.969 ±0.5 10 




0 < 77 <1 


r TBV,pGBV 


0.758 ±0.057 


0.362 ±0.1 05 


0.354±0.101 






^TBV.pGBV 


1.086 ±0.068 


1.455 ± 1.401 


1.251 ±1.310 



Averages and standard errors based on 1 00 replicates of simulated data are listed for prediction accuracy, r pGBV ,TBv. ar) d °i as < ^pgbv.tbv/ of each trait. For the 
settings of prior probability that a SNP has zero effect, tt, see Table 1. 
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trait MCBayes analysis when n =0 while, when n was var- 
ied and inferred, single-trait analysis showed considerably 
higher accuracy than multi-trait MCBayes and varBayes 
analyses as shown in Table 1 and Table 2. The prediction 
with single-trait analysis was almost unbiased In Data I, 
but that was likely to be biased in Data II with &tbv, p gbv 
much deviated from 1 especially when n was varied and 
inferred. 

In Data III with missing phenotypes, which was derived 
from Data I by removing some phenotypic records, r TB v, 
pGbv was decreased for all traits in all methods due to the 
smaller sample size in comparison with Data I (Table 3). 
The rate of decrease in r TB v, P GBV was greater for traits B 
and C than A as the greater reduction of the sample size 
was made in B and C. Multi-trait MCBayes analysis pro- 
vided higher prediction accuracy than multi-trait varBayes 
analysis for all traits except for trait C withzr varied and in- 
ferred. For trait B, multi-trait analysis showed higher pre- 
diction accuracy than single-trait analysis, but gave less 
accurate prediction for trait C. For the bias of predicted 



GBV, &tbv, P GBv» the MCBayes analysis with n =0 was 
almost unbiased for both multi-trait and single-trait 
analyses, where &tbv, p gbv ranged 0.931 to 0.998, although 
the standard error of Z>tbv, p gbv was increased for 
low-heritability traits B and C, for which the MCBayes 
analyses with n varied and varBayes showed much biased 
prediction. 

We listed r TB v, p gbv and 1>tbv, p gbv obtained by cross- 
validation conducted in a population of the 1001st gener- 
ation of Data I as well as r yjpGBV and & y , P GBv> correlation 
between predicted GBV and phenotype and regression of 
phenotype on predicted GBV, in Table 4. The prediction 
accuracy evaluated with cross-validation was considerably 
higher than that evaluated with a test population in the 
next generation, with the prediction bias being similarly 
for both evaluated with cross-validation and use of a test 
population. The relative merits in performance of predic- 
tion between the methods were also similar for both eva- 
luations. It was shown in Table 4 that the relational 
expression, r TB v, P GBV = r YipGBV /h [30], held with h being a 



Table 4 Accuracy and bias of predicted GBVs evaluated with cross-validation in Data I 



Method 




Trait A 


Trait B 


Trait C 


MCBayes 77 = 0 


'TBV.pGBV 


0.832 ±0.039 


0.611 ±0.095 


0.501 ±0.083 




^TBV.pGBV 


1.01 6 ±0.022 


1.072 ±0.231 


1.052 ±0.341 




r y,pGBV 


0.741 ±0.037 


0.191 ±0.045 


0.160 ±0.050 




ky.pGBV 


1.01 3 ±0.01 5 


1.062 ±0.1 25 


1.039 ±0.1 30 


0<77 <1 


'TBV.pGBV 


0.791 ±0.048 


0.603 ±0.1 12 


0.390 ±0.1 27 




^BV.pGBV 


1.1 32 ±0.064 


1.210 ±0.303 


1.1 80 ±0.379 




'y.pGBV 


0.705 ± 0.045 


0.191 ±0.046 


0.121 ±0.060 




by.pGBV 


1.131 ±0.065 


1 .210 ± 0.1 70 


1.1 19 ±0.373 


varBayes 77 = 0 


'TBV.pGBV 


0.81 3 ±0.049 


0.620 ±0.108 


0.470 ±0.1 18 




^TBV.pGBV 


1 .080 ± 0.048 


0.994 ±0.1 57 


0.963 ±0.201 




fy.pGBV 


0.722 ± 0.047 


0.1 87 ±0.056 


0.143 ±0.058 




by.pGBV 


1 .072 ± 0.049 


0.945 ±0.1 95 


0.931 ±0.289 


0 < 77 <1 


'TBV.pGBV 


0.779 ±0.059 


0.593 ±0.1 11 


0.423 ±0.1 23 




^TBV.pGBV 


0.944 ± 0.040 


0.816±0.139 


0.662 ±0.1 53 




r y,pGBV 


0.690 ± 0.056 


0.1 80 ±0.055 


0.1 25 ±0.062 




ky.pGBV 


0.935 ± 0.039 


0.787 ±0.166 


0.626 ± 0.255 


single-trait 77 = 0 


'TBV.pGBV 


0.826 ± 0.040 


0.51 5 ±0.073 


0.505 ± 0.074 


(MCBayes) 


^BV.pGBV 


0.997 ± 0.023 


1.073 ±0.270 


1.030 ±0.303 




'y.pGBV 


0.735 ±0.039 


0.1 59 ±0.045 


0.1 62 ±0.044 




by.pGBV 


0.993 ±0.01 2 


1.071 ±0.162 


1 .055 ± 0.222 


0 < 77 <1 


'TBV.pGBV 


0.821 ±0.039 


0.531 ±0.099 


0.522 ± 0.094 




^TBV.pGBV 


1.131 ±0.046 


1.265 ±0.555 


1.1 92 ±0.405 




'y.pGBV 


0.731 ±0.037 


0.164 ±0.051 


0.1 64 ±0.048 




ky.pGBV 


1.1 27 ±0.043 


1 .249 ± 0.482 


1 .205 ± 0.457 


Averages and standard errors evaluated with 10- fold cross-validation are listed based on 100 replicates of simulated data in Data 1 are listed for prediction 
accuracy, 7 p gbv,tbv an d bias, bpGBv.TBv* as we ll as correlation between phenotypic value and predicted GBV, r yiPGBV , and regression of phenotypic value on 
predicted GBV, b ypG Bv , of each trait. For the settings of prior probability that a SNP has zero effect, n, see Table 1 . 
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square-root of heritability and that ^tbv, p gbv was almost 
the same as & y , P GBV- 

Correlation coefficients between predicted GBVs of trait 
A, B and C were listed in Table 5 for multi-trait MCBayes 
and varBayes analyses and single-trait MCBayes analysis 
in all datasets as well as the correlation coefficients be- 
tween simulated breeding values (TBVs) of three traits in 
the 1002nd generation. Although the correlation of simu- 
lated breeding values between traits A and B was expected 
to be 0.72, the actual correlations obtained were biased 
upwards: 0.755 (0.133) and 0.735 (0.166) in Data I (III) 
and Data II, respectively. Multi-trait analysis could better 
estimate these correlations compared to single-trait ana- 
lysis as seen in the correlation between predicted GBVs 
between traits A and B. For trait C which was genetically 
independent of A and B, the correlations between pre- 
dicted GBVs for traits A and C and traits B and C did not 
significantly deviate from zero (Table 5). 

Discussion 

In this study, we proposed Bayesian methods for simul- 
taneously predicting GBVs for multiple traits, where two 
computational procedures were devised using MCMC it- 
eration and variational approximation, referred to as 
MCBayes and varBayes, respectively. A Bayesian model 



for simultaneously analyzing multiple traits was obtained 
by extending a Bayesian model for single-trait genomic 
selection proposed by [2] and [4] to multi-trait case. We 
introduced the prior probability that a SNP has zero ef- 
fect, n, and accordingly, MCBayes with n fixed at 0, 
meaning that all SNPs are included in a model as covari- 
ates, constructs a model for multi-trait GBV prediction 
in a similar manner to BayesD. On the other hand, the 
MCMC procedure of MCBayes with n variable and in- 
ferred from the data is not exactly the same as that of 
BayesDzr, where a prior distribution of the variance and 
covariance matrix of SNP effects, ~L gb is assumed to be 
independent of whether the SNP is included (yi = 1) or 
excluded(}7 = 0) from the model in MCBayes, as seen in 
(4), while it is dependent on/; and takes different forms 
for yi = 1 and y i = 0 in BayesD7T. This modification allows 
MCMC iteration of MCBayes to be performed with only 
Gibbs sampling. 

In the simultaneous analysis of multiple traits for con- 
structing a GBV prediction model, the computational 
burden greatly increases depending on the number of 
analyzed traits in comparison with single-trait analysis. 
We developed a variational approximation procedure, 
varBayes, for MCBayes to reduce the computational 
time for multi-trait analysis. In varBayes, the joint 



Table 5 Correlations between predicted GBVs for trait A, B and C 



Data Method 


A-B (0.72) 


A-C (0.0) 


B-C (0.0) 


1 MCBayes 77 = 0 


0.588 ±0.1 81 


-0.071 ±0.175 


-0.091 ±0.199 


0<77<1 


0.699 ±0.1 88 


-0.057 ± 0.250 


-0.076 ±0.301 


varBayes 77 = 0 


0.688 ±0.1 84 


-0.1 00 ±0.241 


-0.077 ±0.279 


0<77 <1 


0.644 ±0.1 69 


-0.053 ±0.1 92 


-0.058 ± 0.247 


Single-trait (MCBayes) 77 = 0 


0.446 ±0.1 32 


0.058 ± 0.096 


0.096 ±0.1 21 


0<77 <1 


0.452 ±0.183 


0.035 ± 0.090 


0.060 ±0.1 11 


Simulated BV 


0.755 ±0.1 33 


0.003 ± 0.054 


0.004 ±0.060 


II MCBayes 77 = 0 


0.606 ±0.1 97 


-0.1 29 ±0.1 13 


-0.071 ±0.140 


0<77<1 


0.721 ±0.222 


-0.161 ±0.149 


-0.1 28 ±0.1 79 


varBayes 77 = 0 


0.61 4 ±0.206 


-0.1 76 ±0.1 35 


-0.090 ±0.1 79 


0<77 <1 


0.671 ±0.210 


-0.1 53 ±0.1 68 


-0.047 ±0.1 93 


Single-trait (MCBayes) 77 = 0 


0.394 ±0.1 15 


0.020 ± 0.075 


0.072 ±0.1 11 


0<77<1 


0.479 ± 0.206 


-0.01 8 ±0.084 


0.022 ± 0.094 


Simulated BV 


0.735 ±0.1 66 


-0.031 ±0.054 


-0.01 2 ±0.041 


III MCBayes 77 = 0 


0.530 ±0.1 95 


-0.032 ± 0.223 


-0.037 ±0.238 


0<77 <1 


0.662 ± 0.208 


-0.014 ±0.326 


-0.01 3 ±0.369 


varBayes 77 = 0 


0.555 ±0.1 94 


-0.028 ±0.2 18 


-0.023 ± 0.248 


0<77 <1 


0.455 ±0.1 67 


-0.001 ±0.158 


0.021 ±0.190 


Single-trait (MCBayes) 77 = 0 


0.315 ±0.106 


0.029 ±0.1 05 


0.080 ±0.1 34 


0<77<1 


0.323 ±0.1 38 


0.024 ± 0.097 


0.057 ±0.1 30 



Averages and standard errors are listed based on 100 replicates of simulated data in Data I and Data III and 
simulated breeding values, where expected correlations are 0.72, 0.0 and 0.0 for trait-pairs A-B, A-C and B-C 
between simulated breeding values are the same as those in Data I. 



20 replicates in Data II. Simulated 
as listed in parentheses. In Data II 



BV indicates 
, correlations 
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posterior distribution of parameters was approximated 
by a factorized function, each component of which 
approximated marginal posterior distribution of each 
parameter and was referred to as a variational posterior. 
Variational posteriors were shown to be well-known dis- 
tribution functions such as normal or inverse Wishart 
that could be derived by simple non-MCMC based nu- 
merical iteration. 

In genomic selection, it is important to construct a 
model that enables accurate prediction for GVBs. There- 
fore, precise point estimation of the model parameters is 
more relevant rather than the construction of their poster- 
ior distributions. Accordingly, the evaluation of loss of 
prediction accuracy with varBayes in comparison with 
MCBayes would be suitable for the evaluation of approxi- 
mation accuracy of varBayes. Using simulation experi- 
ments, we investigated the performance of the prediction 
model constructed with multi-trait analysis compared with 
single-trait analysis as well as the model constructed using 
variational approximation. Moreover, the performance of 
multi-trait analysis in the case of missing phenotypic 
records commonly occurring in the treatment of the ac- 
tual data of multiple traits were evaluated based on the 
results of simulations. These points are discussed below 
including the computational time and the possible exten- 
sion of prediction model considering polygenic effects. 

Increase in accuracy for GBV prediction with multi-trait 
analysis 

We evaluated the increase in prediction accuracy with 
multi-trait analysis in comparison with single-trait ana- 
lysis using the datasets without missing phenotypes, 
Data I and Data II (Table 1 and Table 2). For trait A hav- 
ing heritability 0.8, multi-trait analysis and single-trait 
analysis made no difference in the prediction accuracy 
with MCBayes. Therefore, the advantage of multi-trait 
analysis over single-trait analysis in predicting GBV is 
negligible for high-heritability traits. However, we antici- 
pate the increase in accuracy with multi-trait analysis for 
low-heritability traits utilizing correlations with high- 
heritability traits. Actually, for low-heritability trait B 
with heritability 0.1, which has highly correlated with 
trait A (/3gab = 0.72), prediction accuracy was increased 
with multi-trait analysis. As trait C was not genetically 
correlated with trait A, the accuracy of predicting the 
GBV of C was not improved with multi-trait MCBayes 
analysis. The accuracy of predicted GBV of trait B was 
also increased with multi-trait varBayes analysis in Data 
I and Data II in comparison with single-trait MCBayes 
analysis while the prediction accuracy of trait C was 
lower with multi-trait varBayes analysis (Table 1 and 
Table 2). The genetic correlations between traits A and 
B were better estimated with predicted GBVs in multi- 
trait analysis than in single-trait analysis as shown in 



Table 5. Therefore, it can be concluded that low-heritability 
traits (heritability around 0.1) are better predicted for GBVs 
utilizing their correlations with high-heritability traits 
(heritability around 0.8), if any, with multi-trait analysis. 
However, the benefit of multi-trait analysis would be subde 
for high-heritability traits. A similar finding was also 
reported in [13]. For uncorrected low-heritability traits 
such as trait C, it is likely that the prediction using multi- 
trait analysis with n varied and inferred is less effective in 
comparison with single-trait analysis. Therefore, multi-trait 
analysis with jt fixed at 0, which allows highly accurate pre- 
diction for correlated traits while retaining prediction 
accuracy comparable with single-trait analysis for uncorre- 
cted low-heritability traits, would be a suitable method of 
choice for multi-trait genomic selection. 

Approximation accuracy of variational procedure for 
MCMC estimation 

Generally, constructing a GBV prediction model with 
MCMC estimation based on genotypic records for tens 
of thousands of SNPs and phenotypes for hundreds of 
individuals requires considerable computational time 
even for single-trait cases. Much more computational 
burden would be imposed in constructing a model in 
multi-trait analysis, depending on the number of traits 
of interest. Therefore, we proposed a computationally 
cost-effective method, varBayes, approximating MCMC 
based method, MCBayes, using a variational approxima- 
tion procedure. 

Simulation experiments showed that the prediction ac- 
curacy was lower with varBayes than with MCBayes in 
multi-trait analysis but the rate of loss of accuracy was not 
remarkable and was less than 10 percent for traits A and 
B under the same setting of n while it was greater for C 
(Table 1 and Table 2). The prediction accuracy for trait B 
correlated with high-heritability trait A was still higher 
with multi-trait varBayes analysis than with single-trait 
MCBayes analysis indicating that varBayes could well 
utilize the information on the correlation structure in 
multiple traits. Actually, multi-trait varBayes analysis 
could better capture the genetic correlation between A 
and B than single-trait MCbayes analysis (Table 5). 

The computational time was greatly reduced for 
multi-trait varBayes analysis in comparison with multi- 
trait MCBayes analysis. We carried out all computations 
using a Fortran program written to implement multiple- 
trait analysis on a computer having two CPUs each with 
a quad-core processor (Intel Xeon 2.4GHz). In Data I, 
where 100 replicates of datasets each including geno- 
types of 1010 SNPs for 1000 individuals were simulated, 
varBayes took only 12 minutes with n fixed and 22 min- 
utes with 7i varied for 100 times of model constructions 
while the computational time for MCBayes was respect- 
ively 440 and 435 minutes. Therefore, the average 
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computational time required to construct a prediction 
model for each dataset in Data I was less than 15 sec- 
onds for varBayes and was more than 4 minutes for 
MCBayes. For a larger data, Data II, that included 20 
replicates of datasets each consisting of genotypes of 
10100 SNPs for 1000 individuals, the computational 
time was 21 and 19 minutes for varBayes with n fixed 
and ji varied, respectively, and the time for MCBayes 
were increased to more than 12 hours with the average 
computational times for each model construction being 
about 1 minute for varBayes and more than 30 minutes 
for MCBayes. In cross-validation process in which a 
total of one thousand repetitions of model constructions 
were performed in Data I, the computational time was 
about 3days with MCBayes while varBayes completed 
the analysis within only 4 hours. 

Taking computational time and prediction accuracy 
into account, varBayes is considered a useful method for 
multi-trait genomic selection, which can rapidly con- 
struct a prediction model that is less accurate than that 
with the MCMC-based method for multi-trait analysis, 
but is more accurate than that with single-trait analysis 
for correlated traits. The usefulness of varBayes would 
be more remarkable for simultaneous prediction of 
GBVs of a large number of traits based on a huge num- 
ber of SNPs where the application of an MCMC-based 
method might be prohibited. 

Multi-trait analysis of dataset with missing phenotypes 

In Data III, we simulated the datasets under the same 
condition as Data I except that some phenotypes were 
assumed to be unobserved. In short, we assumed that 
phenotypes of traits A, B and C were not available for 
200, 500 and 500 individuals, respectively, in a total of 
1000 individuals with only 100 individuals having the 
phenotypes of all three traits. In multi-trait analysis, 
missing phenotypes of individuals can be estimated with 
their observed phenotypes of other traits using (18), 
which indicates that residual effects of missing pheno- 
types can be restored from those of observed pheno- 
types. When the model fitting is successful for observed 
phenotypes, the residual effects of the phenotypes are 
well estimated by subtracting SNP effects and other 
fixed effects from the phenotypic effects, and those of 
missing phenotypes are suitably obtained by (18) utiliz- 
ing the environmental correlation (covariance) between 
observed and missing phenotypes. Therefore, by assum- 
ing non-zero environmental correlation between traits 
(Peab = 0.1, p EAC = 0.2, pebc = 0.3), the loss of prediction 
accuracy caused by missing phenotypes was anticipated 
to be less with multi-trait analysis than with single-trait 
analysis. We expected that, in Data III, the prediction ac- 
curacy of trait C, environmentally correlated with trait A 
O^eac = 0.2), was maintained higher in the presence of 



missing records for some phenotypes using multi-trait 
analysis in comparison with single-trait analysis. How- 
ever, the rate of loss for the prediction accuracy for trait 
C was similar between multi-trait and single-trait ana- 
lyses as seen in the comparison of the prediction accur- 
acies between Data I and Data III (Tables 1 and 3). The 
utility of implementation of (18) for imputing missing 
phenotypes in the process of multi-trait model construc- 
tion remains unclear in the settings of simulation 
adopted here. 

Model extension by including polygenic effects 

We can modify the Bayesian model (1) by including 
polygenic effects as follows; 

N 

Yi = Xib + J2Yl u ilSi+Vi + ei, (19) 

where v, is a vector of polygenic effects for multiple 
traits and assumed to follow a multivariate normal dis- 
tribution, v, ~7V(0, £„) with E„ being a variance- 
covariance matrix of polygenic effects. The polygenic 
effects for all individuals of a training population and a 
test population are collectively denoted by V, then the 
variance-covariance matrix of V can be expressed as 
A ® £„, where A is additive genetic relationship matrix 
for analyzed individuals, computed from the information 
of pedigree or markers. When the low-density markers 
are used for the analysis, a considerable portion of gen- 
etic effects could not be captured by markers. Accord- 
ingly, if pedigree information is available, inclusion of 
polygenic effects estimated based on pedigree is benefi- 
cial in predicting breeding values for genomic selection. 
The revised model (19) can be similarly treated by 
MCBayes and varBayes as the model (1), where estima- 
tion steps for additional covariates and parameters, V 
and Ev, can be easily implemented in the procedures of 
Bayesian multi-trait model construction using MCMC 
iteration and variational approximation. When no pedi- 
gree information is available, A is computed from mar- 
ker information as a genomic relationship matrix. 
However, the model (19) includes all available markers 
as covariates, resulting in redundantly using the same 
marker information in the model fitting. In the availabil- 
ity of high-density SNP markers, which is the case for 
some species of animals and plants currently, the genetic 
relationships between individuals can be well captured 
by markers themselves in the multi-locus model (1), the 
benefit from the inclusion of the polygenic effects in the 
model would be subtle [31]. 

Conclusion 

In this study, we described a statistical model for Bayesian 
simultaneous prediction of GBVs in genomic selection 
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targeting multiple traits and devised an MCMC-based 
method and a computationally cost-effective method util- 
izing the variational approximation procedure, referred to 
as MCBayes and varBayes, respectively, to estimate para- 
meters included in the model. The results of simulation 
experiments showed that the multi-trait analysis that 
could utilize the correlation structure between traits 
allowed more accurate prediction of GBVs for correlated 
traits compared to single-trait analysis that treated each 
trait separately, where, for low-heritability traits correlated 
with high-heritability traits, the prediction accuracy for 
GBVs was remarkably improved with multi-trait analysis. 
Although the prediction accuracy with varBayes was lower 
than that with MCBayes in multi-trait analysis, the rate of 
loss in accuracy was moderate and the accuracy for corre- 
lated low-heritability traits was still higher with varBayes 
analysis compared to single-trait analysis. Considering the 
benefit of greatly reduced computational time, varBayes 
was considered to be a practical method for predicting 
GBVs in multi-trait genomic selection. 
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